######################################################################################
######################################################################################
#Replication code for:

#Mark T. Buntaine, Ryan Jablonski, Daniel Nielson, Paula Pickering
#SMS Texts on Corruption Help Ugandan Voters Hold Elected Councillors Accountable at the Polls

#Prepared by: Mark Buntaine
#Mark Buntaine contact (as of April 2018): buntaine@bren.ucsb.edu

#Compiled using R Version 3.4.3 (version "Kite-Eating Tree") on Mac running OS X 10.13.4
######################################################################################
######################################################################################

######################################################################################
###Packages
######################################################################################
library(lfe) #Version 2.6-2291
#Note: Contains felm() with demeaning and SE clustering
#Note2: on warnings, see: https://cran.r-project.org/web/packages/lfe/vignettes/identification.pdf
library(stargazer) #Version 5.2.1

######################################################################################
###Functions
######################################################################################
felm.ri <- function(formula, dta, treat.var, rand.ob, rand.ob.info.cols, join.var, sims, ...){
  require(lfe)
  print("1. data and rand.ob must have rows organized in identical order")
  print("2. treat.var should be first right-side entry in formula")
  print("3. join.var must have identical name between data and rand.ob")
  
  ate <- coef(felm(formula, data=dta))[1]
  N <- felm(formula, data=dta)$N
  ate.samp.dist <- rep(NA,sims)
  
  for (i in 1:sims){
    dta[,treat.var] <- rand.ob[rand.ob[,join.var] %in% dta[,join.var], rand.ob.info.cols+i]
    
    ate.samp.dist[i] <- coef(felm(formula, data=dta))[1]
  }
  
  p.two.way <- sum(abs(ate)<abs(ate.samp.dist))/sims
  p.one.way.greater <- sum(ate<ate.samp.dist)/sims
  p.one.way.lesser <- sum(ate>ate.samp.dist)/sims
  se <- sd(ate.samp.dist)
  
  coun <- list("ate" = ate, "ate.samp.dist" = ate.samp.dist, "se"=se, "p.two.way" = p.two.way, "p.one.way.greater" = p.one.way.greater, "p.one.way.lesser" = p.one.way.lesser, "N" = N)
  return(coun)
}

felm.ri2 <- function(formula, dta, treat.var, rand.ob, rand.ob.info.cols, join.var, sims, ...){
  require(lfe)
  print("1. data and rand.ob must have rows organized in identical order")
  print("2. treat.var should be first right-side entry in formula")
  print("3. join.var must have identical name between data and rand.ob")
  
  #Note: have to make this for multiple factor of one randomized treatment
  
  ran.coef.num <- length(unique(dta[,treat.var][!is.na(dta[,treat.var])]))-1 #Gives number of coefficients to keep for crossed treatment indicator
  N <- felm(formula, data=dta)$N
  coef <- coef(felm(formula, data=dta))[1:ran.coef.num]
  coef.samp.dist <- matrix(data=NA, nrow=ran.coef.num, ncol=sims)
  #row.names(coef.samp.dist) <- names(coef)
  
  for (i in 1:sims){
    dta[,treat.var] <- rand.ob[rand.ob[,join.var] %in% dta[,join.var], rand.ob.info.cols+i]
    coef.samp.dist[,i] <- coef(felm(formula, data=dta))[1:ran.coef.num]
  }
  
  se <- apply(coef.samp.dist, 1, sd) #Cannot get SEs off non-randomized parameters
  
  p.two.way <- rep(NA,length(coef))
  p.one.way.greater <- rep(NA,length(coef))
  p.one.way.lesser <- rep(NA,length(coef))
  
  for (i in 1:length(coef)){
    p.two.way[i] <- sum(abs(coef[i])<abs(coef.samp.dist[i,]))/sims
    p.one.way.greater[i] <- sum(coef[i]<coef.samp.dist[i,])/sims
    p.one.way.lesser[i] <- sum(coef[i]>coef.samp.dist[i,])/sims
  }
  
  coun <- list("coef" = coef, "coef.samp.dist" = coef.samp.dist, "se"=se, "p.two.way" = p.two.way, "p.one.way.greater" = p.one.way.greater, "p.one.way.lesser" = p.one.way.lesser, "N" = N)
  return(coun)
}

######################################################################################
###Data input, change file paths as appropriate
######################################################################################

#Update working directory as needed to folder with replication materials
#setwd("D:/google drive/Uganda Vote Choice project/Analysis/Replication/Replication_PNAS/")
setwd("~/Google Drive/Uganda Vote Choice project/Analysis/Replication/Replication_PNAS/")

data <- read.csv("./data/MASTER_Analysis_Anon_180424.csv", stringsAsFactors=FALSE)

Budget_RI <- read.csv("./data/Budget_RI.csv", stringsAsFactors=FALSE)
Budget_RI <- Budget_RI[match(data$id.cleaned, Budget_RI$id.cleaned),] #reordering to match "data"

Density_RI <- read.csv("./data/Density_RI.csv", stringsAsFactors=FALSE)
Density_RI <- Density_RI[match(data$id.cleaned, Density_RI$id.cleaned),] #reordering to match "data"


######################################################################################
###Data setup
######################################################################################
schooling.levels <- c("no_schooling","some_primary_s","completed_prim","some_secondary","completed_seco","some_universit","completed_univ","some_post_grad","completed_mast","refuse_to_answ")
data$r.education <- factor(data$r.How_much_schooling_have_you_co, levels=schooling.levels)

#Remaking "bd.multi.treat2" because of unwanted class conversion in "data" from read.csv()
data$bd.multi.treat2 <- ifelse(data$budget.treat==1 & data$density.treat2==1,"1.1",NA)
data$bd.multi.treat2 <- ifelse(data$budget.treat==1 & data$density.treat2==0,"1.0",data$bd.multi.treat2)
data$bd.multi.treat2 <- ifelse(data$budget.treat==0 & data$density.treat2==1,"0.1",data$bd.multi.treat2)
data$bd.multi.treat2 <- ifelse(data$budget.treat==0 & data$density.treat2==0,"0.0",data$bd.multi.treat2)


######################################################################################
###Budget subsets
######################################################################################

#Prior-defined subgroups
budget.good <- subset(data, budget.actual > budget.prior | (budget.actual==budget.prior & budget.actual>=4))
budget.bad <- subset(data, budget.actual < budget.prior | (budget.actual==budget.prior & budget.actual<=2))

#These are the subgroups not defined by priors, for an extended analysis
#budget.positive <- subset(data, budget.actual>=4)
#budget.negative <- subset(data, budget.actual<=2)

#2. Subset where no individual incumbent switched parties and ran again in 2016 (includes elections w/o incumbent individual)
budget.good.lc5.chair.comp <- subset(budget.good, lc5.chair.competitive==1 & lc5.chair.party.switch==0)
budget.good.lc5.councillor.comp <- subset(budget.good, lc5.councillor.competitive==1 & lc5.councillor.party.switch==0 & is.na(lc5.councillor.redistricted2016))
budget.bad.lc5.chair.comp <- subset(budget.bad, lc5.chair.competitive==1 & lc5.chair.party.switch==0)
budget.bad.lc5.councillor.comp <- subset(budget.bad, lc5.councillor.competitive==1 & lc5.councillor.party.switch==0 & is.na(lc5.councillor.redistricted2016))
#budget.positive.lc5.chair.comp <- subset(budget.positive, lc5.chair.competitive==1 & lc5.chair.party.switch==0)
#budget.positive.lc5.councillor.comp <- subset(budget.positive, lc5.councillor.competitive==1 & lc5.councillor.party.switch==0 & is.na(lc5.councillor.redistricted2016))
#budget.negative.lc5.chair.comp <- subset(budget.negative, lc5.chair.competitive==1 & lc5.chair.party.switch==0)
#budget.negative.lc5.councillor.comp <- subset(budget.negative, lc5.councillor.competitive==1 & lc5.councillor.party.switch==0 & is.na(lc5.councillor.redistricted2016))

#Verified Recipient subgroups
budget.good.lc5.chair.comp.c <- subset(budget.good, lc5.chair.competitive==1 & lc5.chair.party.switch==0 & d.X_11_Over_the_last_several_days=="yes")
budget.good.lc5.councillor.comp.c <- subset(budget.good, lc5.councillor.competitive==1 & lc5.councillor.party.switch==0 & is.na(lc5.councillor.redistricted2016) & d.X_11_Over_the_last_several_days=="yes")
budget.bad.lc5.chair.comp.c <- subset(budget.bad, lc5.chair.competitive==1 & lc5.chair.party.switch==0 & d.X_11_Over_the_last_several_days=="yes")
budget.bad.lc5.councillor.comp.c <- subset(budget.bad, lc5.councillor.competitive==1 & lc5.councillor.party.switch==0 & is.na(lc5.councillor.redistricted2016) & d.X_11_Over_the_last_several_days=="yes")


######################################################################################
###Figure 1: Treatment effects of budget disclosures by news type
######################################################################################

##### Main Effects #####

###Budget / LC5 chair results
chair.good <- felm.ri(lc5.chair.inc.vote ~ budget.treat + lc5.chair.intent | location.id | 0 | lc5.chair.name,
                      dta=budget.good.lc5.chair.comp, treat.var = "budget.treat", rand.ob=Budget_RI, 
                      rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)

chair.bad <- felm.ri(lc5.chair.inc.vote ~ budget.treat + lc5.chair.intent | location.id | 0 | lc5.chair.name, 
                     dta=budget.bad.lc5.chair.comp, treat.var = "budget.treat", rand.ob=Budget_RI, 
                     rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)

chair.good.cov <- felm.ri(lc5.chair.inc.vote ~ budget.treat + lc5.chair.intent + b.Q1_living_conditions + b.Q21b_Trust_in_Twaweza_for_info + b.Q24_Powerful_ppl_learning_how_ + b.Q25_Will_counting_votes_be_fai + r.What_is_your_gender + r.education + r.What_is_your_age | location.id | 0 | lc5.chair.name,
                          dta=budget.good.lc5.chair.comp, treat.var = "budget.treat", rand.ob=Budget_RI, 
                          rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)

chair.bad.cov <- felm.ri(lc5.chair.inc.vote ~ budget.treat + lc5.chair.intent + b.Q1_living_conditions + b.Q21b_Trust_in_Twaweza_for_info + b.Q24_Powerful_ppl_learning_how_ + b.Q25_Will_counting_votes_be_fai + r.What_is_your_gender + r.education + r.What_is_your_age | location.id | 0 | lc5.chair.name, 
                         dta=budget.bad.lc5.chair.comp, treat.var = "budget.treat", rand.ob=Budget_RI, 
                         rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)


###Budget / LC5 councillor results
coun.good <- felm.ri(lc5.councillor.inc.vote ~ budget.treat + lc5.councillor.intent | location.id | 0 | lc5.councillor.name,
                     dta=budget.good.lc5.councillor.comp, treat.var = "budget.treat", rand.ob=Budget_RI, 
                     rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)

coun.bad <- felm.ri(lc5.councillor.inc.vote ~ budget.treat + lc5.councillor.intent | location.id | 0 | lc5.councillor.name, 
                    dta=budget.bad.lc5.councillor.comp, treat.var = "budget.treat", rand.ob=Budget_RI, 
                    rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)

coun.good.cov <- felm.ri(lc5.councillor.inc.vote ~ budget.treat + lc5.councillor.intent + b.Q1_living_conditions + b.Q21b_Trust_in_Twaweza_for_info + b.Q24_Powerful_ppl_learning_how_ + b.Q25_Will_counting_votes_be_fai + r.What_is_your_gender + r.education + r.What_is_your_age | location.id | 0 | lc5.councillor.name,
                         dta=budget.good.lc5.councillor.comp, treat.var = "budget.treat", rand.ob=Budget_RI, 
                         rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)

coun.bad.cov <- felm.ri(lc5.councillor.inc.vote ~ budget.treat + lc5.councillor.intent + b.Q1_living_conditions + b.Q21b_Trust_in_Twaweza_for_info + b.Q24_Powerful_ppl_learning_how_ + b.Q25_Will_counting_votes_be_fai + r.What_is_your_gender + r.education + r.What_is_your_age | location.id | 0 | lc5.councillor.name, 
                        dta=budget.bad.lc5.councillor.comp, treat.var = "budget.treat", rand.ob=Budget_RI, 
                        rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)

##### Verified Recipient Effects #####

###Budget / LC5 chair results, verified compliers
chair.good.c <- felm.ri(lc5.chair.inc.vote ~ budget.treat + lc5.chair.intent | location.id | 0 | lc5.chair.name,
                        dta=budget.good.lc5.chair.comp.c, treat.var = "budget.treat", rand.ob=Budget_RI, 
                        rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)

chair.bad.c <- felm.ri(lc5.chair.inc.vote ~ budget.treat + lc5.chair.intent | location.id | 0 | lc5.chair.name, 
                       dta=budget.bad.lc5.chair.comp.c, treat.var = "budget.treat", rand.ob=Budget_RI, 
                       rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)

chair.good.cov.c <- felm.ri(lc5.chair.inc.vote ~ budget.treat + lc5.chair.intent + b.Q1_living_conditions + b.Q21b_Trust_in_Twaweza_for_info + b.Q24_Powerful_ppl_learning_how_ + b.Q25_Will_counting_votes_be_fai + r.What_is_your_gender + r.education + r.What_is_your_age | location.id | 0 | lc5.chair.name,
                            dta=budget.good.lc5.chair.comp.c, treat.var = "budget.treat", rand.ob=Budget_RI, 
                            rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)

chair.bad.cov.c <- felm.ri(lc5.chair.inc.vote ~ budget.treat + lc5.chair.intent + b.Q1_living_conditions + b.Q21b_Trust_in_Twaweza_for_info + b.Q24_Powerful_ppl_learning_how_ + b.Q25_Will_counting_votes_be_fai + r.What_is_your_gender + r.education + r.What_is_your_age | location.id | 0 | lc5.chair.name, 
                           dta=budget.bad.lc5.chair.comp.c, treat.var = "budget.treat", rand.ob=Budget_RI, 
                           rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)


###Budget / LC5 councillor results, verified compliers
coun.good.c <- felm.ri(lc5.councillor.inc.vote ~ budget.treat + lc5.councillor.intent | location.id | 0 | lc5.councillor.name,
                       dta=budget.good.lc5.councillor.comp.c, treat.var = "budget.treat", rand.ob=Budget_RI, 
                       rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)

coun.bad.c <- felm.ri(lc5.councillor.inc.vote ~ budget.treat + lc5.councillor.intent | location.id | 0 | lc5.councillor.name, 
                      dta=budget.bad.lc5.councillor.comp.c, treat.var = "budget.treat", rand.ob=Budget_RI, 
                      rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)

coun.good.cov.c <- felm.ri(lc5.councillor.inc.vote ~ budget.treat + lc5.councillor.intent + b.Q1_living_conditions + b.Q21b_Trust_in_Twaweza_for_info + b.Q24_Powerful_ppl_learning_how_ + b.Q25_Will_counting_votes_be_fai + r.What_is_your_gender + r.education + r.What_is_your_age | location.id | 0 | lc5.councillor.name,
                           dta=budget.good.lc5.councillor.comp.c, treat.var = "budget.treat", rand.ob=Budget_RI, 
                           rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)

coun.bad.cov.c <- felm.ri(lc5.councillor.inc.vote ~ budget.treat + lc5.councillor.intent + b.Q1_living_conditions + b.Q21b_Trust_in_Twaweza_for_info + b.Q24_Powerful_ppl_learning_how_ + b.Q25_Will_counting_votes_be_fai + r.What_is_your_gender + r.education + r.What_is_your_age | location.id | 0 | lc5.councillor.name, 
                          dta=budget.bad.lc5.councillor.comp.c, treat.var = "budget.treat", rand.ob=Budget_RI, 
                          rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)


###Main Effects, Piped Values
te.good <- c(chair.good$ate, chair.good.cov$ate, coun.good$ate, coun.good.cov$ate)
se.good <- c(chair.good$se, chair.good.cov$se, coun.good$se, coun.good.cov$se)
base.c.good <- c(mean(budget.good.lc5.chair.comp$lc5.chair.inc.vote[budget.good.lc5.chair.comp$budget.treat==0], na.rm=TRUE),
                 mean(budget.good.lc5.chair.comp$lc5.chair.inc.vote[budget.good.lc5.chair.comp$budget.treat==0], na.rm=TRUE),
                 mean(budget.good.lc5.councillor.comp$lc5.councillor.inc.vote[budget.good.lc5.councillor.comp$budget.treat==0], na.rm=TRUE),
                 mean(budget.good.lc5.councillor.comp$lc5.councillor.inc.vote[budget.good.lc5.councillor.comp$budget.treat==0], na.rm=TRUE))
n.good <- c(chair.good$N,coun.good$N)

te.bad <- c(chair.bad$ate, chair.bad.cov$ate, coun.bad$ate, coun.bad.cov$ate)
se.bad <- c(chair.bad$se, chair.bad.cov$se, coun.bad$se, coun.bad.cov$se)
base.c.bad <- c(mean(budget.bad.lc5.chair.comp$lc5.chair.inc.vote[budget.bad.lc5.chair.comp$budget.treat==0], na.rm=TRUE),
                mean(budget.bad.lc5.chair.comp$lc5.chair.inc.vote[budget.bad.lc5.chair.comp$budget.treat==0], na.rm=TRUE),
                mean(budget.bad.lc5.councillor.comp$lc5.councillor.inc.vote[budget.bad.lc5.councillor.comp$budget.treat==0], na.rm=TRUE),
                mean(budget.bad.lc5.councillor.comp$lc5.councillor.inc.vote[budget.bad.lc5.councillor.comp$budget.treat==0], na.rm=TRUE))
n.bad <- c(chair.bad$N,coun.bad$N)

###Compliers, Piped Valued
te.good.c <- c(chair.good.c$ate, chair.good.cov.c$ate, coun.good.c$ate, coun.good.cov.c$ate)
se.good.c <- c(chair.good.c$se, chair.good.cov.c$se, coun.good.c$se, coun.good.cov.c$se)
base.c.good.c <- c(mean(budget.good.lc5.chair.comp.c$lc5.chair.inc.vote[budget.good.lc5.chair.comp.c$budget.treat==0], na.rm=TRUE),
                   mean(budget.good.lc5.chair.comp.c$lc5.chair.inc.vote[budget.good.lc5.chair.comp.c$budget.treat==0], na.rm=TRUE),
                   mean(budget.good.lc5.councillor.comp.c$lc5.councillor.inc.vote[budget.good.lc5.councillor.comp.c$budget.treat==0], na.rm=TRUE),
                   mean(budget.good.lc5.councillor.comp.c$lc5.councillor.inc.vote[budget.good.lc5.councillor.comp.c$budget.treat==0], na.rm=TRUE))
n.good.c <- c(chair.good.c$N,coun.good.c$N)

te.bad.c <- c(chair.bad.c$ate, chair.bad.cov.c$ate, coun.bad.c$ate, coun.bad.cov.c$ate)
se.bad.c <- c(chair.bad.c$se, chair.bad.cov.c$se, coun.bad.c$se, coun.bad.cov.c$se)
base.c.bad.c <- c(mean(budget.bad.lc5.chair.comp.c$lc5.chair.inc.vote[budget.bad.lc5.chair.comp.c$budget.treat==0], na.rm=TRUE),
                  mean(budget.bad.lc5.chair.comp.c$lc5.chair.inc.vote[budget.bad.lc5.chair.comp.c$budget.treat==0], na.rm=TRUE),
                  mean(budget.bad.lc5.councillor.comp.c$lc5.councillor.inc.vote[budget.bad.lc5.councillor.comp.c$budget.treat==0], na.rm=TRUE),
                  mean(budget.bad.lc5.councillor.comp.c$lc5.councillor.inc.vote[budget.bad.lc5.councillor.comp.c$budget.treat==0], na.rm=TRUE))
n.bad.c <- c(chair.bad.c$N,coun.bad.c$N)

##### Combined Plot of Main / Verified Recipient Effects  #####

#pdf("VoteChoice_Main_171122.pdf",height=7, width=6)
heights <- c(8,6,3,1)
layout(matrix(1:5, ncol = 1), widths = 1, heights = c(4,4,4,4,1), respect = FALSE)

###Good News Plot (Main)
par(mar=c(3.1,9.1,3.1,6.1), mgp=c(2,0.8,0), cex.main=1.6, cex.axis=1.4, bty="o")
plot(1, type="n", xlim=c(-0.1,0.1), ylim=c(0,9), yaxt="n", ylab="", xlab="Treatment Effect", cex.lab=1.4)
title("Good News (All Subjects)", line=0.5)
points(x=te.good,y=heights, pch=c(1,16,1,16), cex=3)
abline(v=0, lty=2)

for (i in 1:length(te.good)){
  lines(x=c(te.good[i]-1.96*se.good[i],te.good[i]+1.96*se.good[i]), y=c(heights[i],heights[i]), lwd=1)
  lines(x=c(te.good[i]-1.645*se.good[i],te.good[i]+1.645*se.good[i]), y=c(heights[i],heights[i]), lwd=3)
}

axis(2, at=c(7,2), labels=c("LC5 Chair", "LC5 Councillor"), las=2, tck=0, cex.axis=1.4)

axis(4, at=c(8.1,3.1), labels=c(paste("n=",n.good[1],sep=""), paste("n=",n.good[2],sep="")), las=2, tck=0, cex.axis=1.4)
axis(4, at=6.2, labels=substitute(bar(y[c])==b1,list(b1=round(base.c.good[1],2))),
     las=2, tck=0, cex.axis=1.4)
axis(4, at=1.2, labels=substitute(bar(y[c])==b3,list(b3=round(base.c.good[3],2))),
     las=2, tck=0, cex.axis=1.4)

###Good News Plot (Verified Treated)
par(mar=c(3.1,9.1,3.1,6.1), mgp=c(2,0.8,0), cex.main=1.6, cex.axis=1.4, bty="o")
plot(1, type="n", xlim=c(-0.1,0.1), ylim=c(0,9), yaxt="n", ylab="", xlab="Treatment Effect", cex.lab=1.4)
title("Good News (Verified Recipients)", line=0.5)
points(x=te.good.c,y=heights, pch=c(1,16,1,16), cex=3)
abline(v=0, lty=2)

for (i in 1:length(te.good)){
  lines(x=c(te.good.c[i]-1.96*se.good.c[i],te.good.c[i]+1.96*se.good.c[i]), y=c(heights[i],heights[i]), lwd=1)
  lines(x=c(te.good.c[i]-1.645*se.good.c[i],te.good.c[i]+1.645*se.good.c[i]), y=c(heights[i],heights[i]), lwd=3)
}

axis(2, at=c(7,2), labels=c("LC5 Chair", "LC5 Councillor"), las=2, tck=0, cex.axis=1.4)

axis(4, at=c(8.1,3.1), labels=c(paste("n=",n.good.c[1],sep=""), paste("n=",n.good.c[2],sep="")), las=2, tck=0, cex.axis=1.4)
axis(4, at=6.2, labels=substitute(bar(y[c])==b1,list(b1=round(base.c.good.c[1],2))),
     las=2, tck=0, cex.axis=1.4)
axis(4, at=1.2, labels=substitute(bar(y[c])==b3,list(b3=round(base.c.good.c[3],2))),
     las=2, tck=0, cex.axis=1.4)


###Bad News Plot (Main)
par(mar=c(3.1,9.1,3.1,6.1), mgp=c(2,0.8,0), cex.main=1.6, cex.axis=1.4, bty="o")
plot(1, type="n", xlim=c(-0.1,0.1), ylim=c(0,9), yaxt="n", ylab="", xlab="Treatment Effect", cex.lab=1.4)
title("Bad News (All Subjects)", line=0.5)
points(x=te.bad,y=heights, pch=c(1,16,1,16), cex=3)
abline(v=0, lty=2)


for (i in 1:length(te.bad)){
  lines(x=c(te.bad[i]-1.96*se.bad[i],te.bad[i]+1.96*se.bad[i]), y=c(heights[i],heights[i]), lwd=1)
  lines(x=c(te.bad[i]-1.645*se.bad[i],te.bad[i]+1.645*se.bad[i]), y=c(heights[i],heights[i]), lwd=3)
}

axis(2, at=c(7,2), labels=c("LC5 Chair", "LC5 Councillor"), las=2, tck=0, cex.axis=1.4)

axis(4, at=c(8.1,3.1), labels=c(paste("n=",n.bad[1],sep=""), paste("n=",n.bad[2],sep="")), las=2, tck=0, cex.axis=1.4)
axis(4, at=6.2, labels=substitute(bar(y[c])==b1,list(b1=round(base.c.bad[1],2))),
     las=2, tck=0, cex.axis=1.4)
axis(4, at=1.2, labels=substitute(bar(y[c])==b3,list(b3=round(base.c.bad[3],2))),
     las=2, tck=0, cex.axis=1.4)

###Bad News Plot (Verified Treated)
par(mar=c(3.1,9.1,3.1,6.1), mgp=c(2,0.8,0), cex.main=1.6, cex.axis=1.4, bty="o")
plot(1, type="n", xlim=c(-0.1,0.1), ylim=c(0,9), yaxt="n", ylab="", xlab="Treatment Effect", cex.lab=1.4)
title("Bad News (Verified Recipients)", line=0.5)
points(x=te.bad.c,y=heights, pch=c(1,16,1,16), cex=3)
abline(v=0, lty=2)

for (i in 1:length(te.bad)){
  lines(x=c(te.bad.c[i]-1.96*se.bad.c[i],te.bad.c[i]+1.96*se.bad.c[i]), y=c(heights[i],heights[i]), lwd=1)
  lines(x=c(te.bad.c[i]-1.645*se.bad.c[i],te.bad.c[i]+1.645*se.bad.c[i]), y=c(heights[i],heights[i]), lwd=3)
}

axis(2, at=c(7,2), labels=c("LC5 Chair", "LC5 Councillor"), las=2, tck=0, cex.axis=1.4)

axis(4, at=c(8.1,3.1), labels=c(paste("n=",n.bad.c[1],sep=""), paste("n=",n.bad.c[2],sep="")), las=2, tck=0, cex.axis=1.4)
axis(4, at=6.2, labels=substitute(bar(y[c])==b1,list(b1=round(base.c.bad.c[1],2))),
     las=2, tck=0, cex.axis=1.4)
axis(4, at=1.2, labels=substitute(bar(y[c])==b3,list(b3=round(base.c.bad.c[3],2))),
     las=2, tck=0, cex.axis=1.4)

###Legend
par(mar=c(0.1,0.1,0.1,0.1), bty="n")
plot(1, type="n", xlim=c(0,10), ylim=c(0,2), xaxt="n", yaxt="n", ylab="", xlab="")
points(x=c(2.6,5.6),y=c(0.8,0.8), pch=c(1,16), cex=3)
text(x=2.7, y=0.8, pos=4, "Without Covariates", cex=1.4)
text(x=5.7, y=0.8, pos=4, "With Covariates", cex=1.4)

#dev.off()


######################################################################################
###Figure 2: The treatment effect of higher treatment density among treated subjects
######################################################################################

h11a.b.1 <- felm.ri(lc5.chair.inc.vote ~ density.treat2 + lc5.chair.intent + b.Q1_living_conditions + b.Q21b_Trust_in_Twaweza_for_info + b.Q24_Powerful_ppl_learning_how_ + b.Q25_Will_counting_votes_be_fai + r.What_is_your_gender + r.education + r.What_is_your_age | d.block | 0 | lc5.chair.name, 
                    dta=budget.good.lc5.chair.comp[budget.good.lc5.chair.comp$budget.treat==1,],
                    treat.var = "density.treat2", rand.ob=Density_RI,
                    rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)

h11a.b.2 <- felm.ri(lc5.councillor.inc.vote ~ density.treat2 + lc5.councillor.intent + b.Q1_living_conditions + b.Q21b_Trust_in_Twaweza_for_info + b.Q24_Powerful_ppl_learning_how_ + b.Q25_Will_counting_votes_be_fai + r.What_is_your_gender + r.education + r.What_is_your_age | d.block | 0 | lc5.councillor.name, 
                    dta=budget.good.lc5.councillor.comp[budget.good.lc5.councillor.comp$budget.treat==1,],
                    treat.var = "density.treat2", rand.ob=Density_RI,
                    rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)

h11b.b.1 <- felm.ri(lc5.chair.inc.vote ~ density.treat2 + lc5.chair.intent + b.Q1_living_conditions + b.Q21b_Trust_in_Twaweza_for_info + b.Q24_Powerful_ppl_learning_how_ + b.Q25_Will_counting_votes_be_fai + r.What_is_your_gender + r.education + r.What_is_your_age | d.block | 0 | lc5.chair.name, 
                    dta=budget.bad.lc5.chair.comp[budget.bad.lc5.chair.comp$budget.treat==1,],
                    treat.var = "density.treat2", rand.ob=Density_RI,
                    rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)

h11b.b.2 <- felm.ri(lc5.councillor.inc.vote ~ density.treat2 + lc5.councillor.intent + b.Q1_living_conditions + b.Q21b_Trust_in_Twaweza_for_info + b.Q24_Powerful_ppl_learning_how_ + b.Q25_Will_counting_votes_be_fai + r.What_is_your_gender + r.education + r.What_is_your_age | d.block | 0 | lc5.councillor.name, 
                    dta=budget.bad.lc5.councillor.comp[budget.bad.lc5.councillor.comp$budget.treat==1,],
                    treat.var = "density.treat2", rand.ob=Density_RI,
                    rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)

###Piped Values for Density Effect Among Treated, Figure 5
te.good <- c(h11a.b.1$ate, NA, h11a.b.2$ate, NA) #-0.03465979,NA,-0.04034953,NA
se.good <- c(h11a.b.1$se, NA, h11a.b.2$se, NA) #0.02791619,NA,0.05255581,NA
base.c.good <- c(mean(budget.good.lc5.chair.comp$lc5.chair.inc.vote[budget.good.lc5.chair.comp$budget.treat==1 & budget.good.lc5.chair.comp$density.treat2==0], na.rm=TRUE),
                 mean(budget.good.lc5.chair.comp$lc5.chair.inc.vote[budget.good.lc5.chair.comp$budget.treat==1 & budget.good.lc5.chair.comp$density.treat2==0], na.rm=TRUE),
                 mean(budget.good.lc5.councillor.comp$lc5.councillor.inc.vote[budget.good.lc5.councillor.comp$budget.treat==1 & budget.good.lc5.councillor.comp$density.treat==0], na.rm=TRUE),
                 mean(budget.good.lc5.councillor.comp$lc5.councillor.inc.vote[budget.good.lc5.councillor.comp$budget.treat==1 & budget.good.lc5.councillor.comp$density.treat==0], na.rm=TRUE))
n.good <- c(1703,1335)

te.bad <- c(h11b.b.1$ate, NA, h11b.b.2$ate, NA) #0.001373277,NA,-0.035120041,NA
se.bad <- c(h11b.b.1$se, NA, h11b.b.2$se, NA) #0.02863102,NA,0.05652717,NA
base.c.bad <- c(mean(budget.bad.lc5.chair.comp$lc5.chair.inc.vote[budget.bad.lc5.chair.comp$budget.treat==1 & budget.bad.lc5.chair.comp$density.treat2==0], na.rm=TRUE),
                mean(budget.bad.lc5.chair.comp$lc5.chair.inc.vote[budget.bad.lc5.chair.comp$budget.treat==1 & budget.bad.lc5.chair.comp$density.treat2==0], na.rm=TRUE),
                mean(budget.bad.lc5.councillor.comp$lc5.councillor.inc.vote[budget.bad.lc5.councillor.comp$budget.treat==1 & budget.bad.lc5.councillor.comp$density.treat2==0], na.rm=TRUE),
                mean(budget.bad.lc5.councillor.comp$lc5.councillor.inc.vote[budget.bad.lc5.councillor.comp$budget.treat==1 & budget.bad.lc5.councillor.comp$density.treat2==0], na.rm=TRUE))
n.bad <- c(1325,1157)
x.axis.lab <- "Proportion Voting for Incumbent or Incumbent Party"

#pdf("VoteChoice_Fig2_Density_171127.pdf", height=3.7, width=6)
heights <- c(8,6,3,1)
layout(matrix(1:3, ncol = 1), widths = 1, heights = c(4,4,1), respect = FALSE)

###Good News Plot
par(mar=c(3.1,9.1,3.1,6.1), mgp=c(2,0.8,0), cex.main=1.6, cex.axis=1.4, bty="o")
plot(1, type="n", xlim=c(-0.15,0.1), ylim=c(0,9), yaxt="n", ylab="", xlab="Treatment Effect", cex.lab=1.4)
title("Good News", line=0.5)
points(x=te.good,y=heights-1, pch=c(16,1,16,1), cex=3)
abline(v=0, lty=2)

for (i in 1:length(te.good)){
  lines(x=c(te.good[i]-1.96*se.good[i],te.good[i]+1.96*se.good[i]), y=c(heights[i]-1,heights[i]-1), lwd=1)
  lines(x=c(te.good[i]-1.645*se.good[i],te.good[i]+1.645*se.good[i]), y=c(heights[i]-1,heights[i]-1), lwd=3)
}

axis(2, at=c(7,2), labels=c("LC5 Chair", "LC5 Councillor"), las=2, tck=0, cex.axis=1.4)

axis(4, at=c(8.1,3.1), labels=c(paste("n=",n.good[1],sep=""), paste("n=",n.good[2],sep="")), las=2, tck=0, cex.axis=1.4)
axis(4, at=6.2, labels=substitute(bar(y[L])==b1,list(b1=round(base.c.good[1],2))),
     las=2, tck=0, cex.axis=1.4)
axis(4, at=1.2, labels=substitute(bar(y[L])==b3,list(b3=round(base.c.good[3],2))),
     las=2, tck=0, cex.axis=1.4)

###Bad News Plot
par(mar=c(3.1,9.1,3.1,6.1), mgp=c(2,0.8,0), cex.main=1.6, cex.axis=1.4, bty="o")
plot(1, type="n", xlim=c(-0.15,0.1), ylim=c(0,9), yaxt="n", ylab="", xlab="Treatment Effect", cex.lab=1.4)
title("Bad News", line=0.5)
points(x=te.bad,y=heights-1, pch=c(16,1,16,1), cex=3)
abline(v=0, lty=2)

for (i in 1:length(te.bad)){
  lines(x=c(te.bad[i]-1.96*se.bad[i],te.bad[i]+1.96*se.bad[i]), y=c(heights[i]-1,heights[i]-1), lwd=1)
  lines(x=c(te.bad[i]-1.645*se.bad[i],te.bad[i]+1.645*se.bad[i]), y=c(heights[i]-1,heights[i]-1), lwd=3)
}

axis(2, at=c(7,2), labels=c("LC5 Chair", "LC5 Councillor"), las=2, tck=0, cex.axis=1.4)

axis(4, at=c(8.1,3.1), labels=c(paste("n=",n.bad[1],sep=""), paste("n=",n.bad[2],sep="")), las=2, tck=0, cex.axis=1.4)
axis(4, at=6.2, labels=substitute(bar(y[L])==b1,list(b1=round(base.c.bad[1],2))),
     las=2, tck=0, cex.axis=1.4)
#axis(4, at=1.2, labels=substitute(bar(y[L])==b3,list(b3=round(base.c.bad[3],2))),
axis(4, at=1.2, labels=substitute(bar(y[L])==b3,list(b3="0.70")),
     las=2, tck=0, cex.axis=1.4)

###Legend
par(mar=c(0.1,0.1,0.1,0.1), bty="n")
plot(1, type="n", xlim=c(0,10), ylim=c(0,2), xaxt="n", yaxt="n", ylab="", xlab="")
points(x=c(4.5),y=c(0.88), pch=16, cex=3)
text(x=4.6, y=0.8, pos=4, "With Covariates", cex=1.4)

#dev.off()


######################################################################################
###Table 2: Effects of budget treatment and treatment density
######################################################################################

b.ri <- data.matrix(Budget_RI)
d.ri <- data.matrix(Density_RI)
#Note: NAs introduced by coercion warning generated because "district" and "village" variables cannot be made numeric
#Note: these NAs are excluded when only the treatment variables are pasted into the final "BD_RI" dataframe in the block below

BD_RI <- matrix(paste(b.ri,d.ri,sep="."),nrow=16083,ncol=10003)
BD_RI[BD_RI=="0.NA" | BD_RI=="1.NA"] <- NA
BD_RI[BD_RI=="0. 0"] <- "0.0"
BD_RI[BD_RI=="0. 1"] <- "0.1"
BD_RI[BD_RI=="1. 0"] <- "1.0"
BD_RI[BD_RI=="1. 1"] <- "1.1"
BD_RI <- as.data.frame(BD_RI)
BD_RI <- cbind(Budget_RI[,1:3],BD_RI[,4:ncol(BD_RI)])
remove(b.ri,d.ri)

h11a.b.1 <- felm(lc5.chair.inc.vote ~ bd.multi.treat + lc5.chair.intent + b.Q1_living_conditions + b.Q21b_Trust_in_Twaweza_for_info + b.Q24_Powerful_ppl_learning_how_ + b.Q25_Will_counting_votes_be_fai + r.What_is_your_gender + r.education + r.What_is_your_age | d.block | 0 | lc5.chair.name, data=budget.good.lc5.chair.comp)
h11a.b.2 <- felm(lc5.councillor.inc.vote ~ bd.multi.treat + lc5.councillor.intent + b.Q1_living_conditions + b.Q21b_Trust_in_Twaweza_for_info + b.Q24_Powerful_ppl_learning_how_ + b.Q25_Will_counting_votes_be_fai + r.What_is_your_gender + r.education + r.What_is_your_age | d.block | 0 | lc5.councillor.name, data=budget.good.lc5.councillor.comp)
h11b.b.1 <- felm(lc5.chair.inc.vote ~ bd.multi.treat + lc5.chair.intent + b.Q1_living_conditions + b.Q21b_Trust_in_Twaweza_for_info + b.Q24_Powerful_ppl_learning_how_ + b.Q25_Will_counting_votes_be_fai + r.What_is_your_gender + r.education + r.What_is_your_age | d.block | 0 | lc5.chair.name, data=budget.bad.lc5.chair.comp)
h11b.b.2 <- felm(lc5.councillor.inc.vote ~ bd.multi.treat + lc5.councillor.intent + b.Q1_living_conditions + b.Q21b_Trust_in_Twaweza_for_info + b.Q24_Powerful_ppl_learning_how_ + b.Q25_Will_counting_votes_be_fai + r.What_is_your_gender + r.education + r.What_is_your_age | d.block | 0 | lc5.councillor.name, data=budget.bad.lc5.councillor.comp)
#Note: treatment coefficients displayed in a different order than output table

chair.good.crossed <- felm.ri2(lc5.chair.inc.vote ~ bd.multi.treat2 + lc5.chair.intent + b.Q1_living_conditions + b.Q21b_Trust_in_Twaweza_for_info + b.Q24_Powerful_ppl_learning_how_ + b.Q25_Will_counting_votes_be_fai + r.What_is_your_gender + r.education + r.What_is_your_age | d.block | 0 | lc5.chair.name,
                               dta=budget.good.lc5.chair.comp, treat.var = "bd.multi.treat2", rand.ob=BD_RI,
                               rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)


coun.good.crossed <- felm.ri2(lc5.councillor.inc.vote ~ bd.multi.treat2 + lc5.councillor.intent + b.Q1_living_conditions + b.Q21b_Trust_in_Twaweza_for_info + b.Q24_Powerful_ppl_learning_how_ + b.Q25_Will_counting_votes_be_fai + r.What_is_your_gender + r.education + r.What_is_your_age | d.block | 0 | lc5.councillor.name, 
                              dta=budget.good.lc5.councillor.comp, treat.var = "bd.multi.treat2", rand.ob=BD_RI,
                              rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)

chair.bad.crossed <- felm.ri2(lc5.chair.inc.vote ~ bd.multi.treat2 + lc5.chair.intent + b.Q1_living_conditions + b.Q21b_Trust_in_Twaweza_for_info + b.Q24_Powerful_ppl_learning_how_ + b.Q25_Will_counting_votes_be_fai + r.What_is_your_gender + r.education + r.What_is_your_age | d.block | 0 | lc5.chair.name,
                              dta=budget.bad.lc5.chair.comp, treat.var = "bd.multi.treat2", rand.ob=BD_RI,
                              rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)


coun.bad.crossed <- felm.ri2(lc5.councillor.inc.vote ~ bd.multi.treat2 + lc5.councillor.intent + b.Q1_living_conditions + b.Q21b_Trust_in_Twaweza_for_info + b.Q24_Powerful_ppl_learning_how_ + b.Q25_Will_counting_votes_be_fai + r.What_is_your_gender + r.education + r.What_is_your_age | d.block | 0 | lc5.councillor.name, 
                             dta=budget.bad.lc5.councillor.comp, treat.var = "bd.multi.treat2", rand.ob=BD_RI,
                             rand.ob.info.cols = 4, join.var = "id.cleaned", sims=9999)

stargazer(h11a.b.1, h11a.b.2, h11b.b.1, h11b.b.2, type = "latex", 
          keep = c("bd.multi.treatbudget0.density1","bd.multi.treatbudget1.density0","bd.multi.treatbudget1.density1"),
          dep.var.labels = c("LC5 Chair", "LC5 Councillor", "LC5 Chair", "LC5 Councillor"),
          column.labels = c("Good News","Bad News"), column.separate = c(2,2),
          dep.var.caption  = "DV: Vote Choice for the Incumbent",
          covariate.labels = c("Treated, Low Density (RI)", "Control, High Density (RI)", "Treated, High Density (RI)"),
          notes = "Note: SEs derived from RI; one-tailed tests as pre-registered; contested elections only",
          df = FALSE, omit.stat = c("rsq","ser"),
          se = list(chair.good.crossed$se,coun.good.crossed$se,chair.bad.crossed$se,coun.bad.crossed$se),
          p = list(chair.good.crossed$p.one.way.greater, coun.good.crossed$p.one.way.greater, chair.bad.crossed$p.one.way.lesser, coun.bad.crossed$p.one.way.lesser),
          report = c('vcsp'), order=c(2,1,3),
          notes.append = FALSE,
          notes.label = "",
          add.lines = list(c("Paired village fixed effects", "Yes", "Yes","Yes","Yes"),
                           c("Covariates", "Yes", "Yes","Yes","Yes"))
)
#Note: stargazer works best with model objects, providing them here to minimize the complexity of the code

